Computing the multifractal spectrum from time series: An algorithmic approach 
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We show that the existing methods for computing the /(«) spectrum from a time series can be 
improved by using a new algorithmic scheme. The scheme relies on the basic idea that the smooth 
convex profile of a typical f(a) spectrum can be fitted with an analytic function involving a set of 
four independant parameters. While the standard existing schemes pH , l20j generally compute only 
an incomplete f(a) spectrum (usually the top portion), we show that this can be overcome by an 
algorithmic approach which is automated to compute the D q and f(a) spectrum from a time series 
for any embedding dimension. The scheme is first tested with the logistic attractor with known 
f(a) curve and subsequently applied to higher dimensional cases. We also show that the scheme 
can be effectively adapted for analysing practcal time series involving noise, with examples from 
two widely different real world systems. Moreover, some preliminary results indicating that the set 
of four independant parameters may be used as diagnostic measures is also included. 

PACS numbers: 05.45.Ac, 05.45.Tp, 05.45.Df 
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It is now well established that multifractal sets 
and objects abound in Nature. A characteristic 
feature of these objects is the self similarity since 
their formation is governed by subtle scaling laws. 
An important tool to analyse these sets is the f(a) 
spectrum which describes how the fractal dimen- 
sions of the interwoven sets with defnite singular- 
ity strength are distributed. In the recent issue 
of Chaos, participating in the discussion "Is the 
normal heart rate chaotic?", many authors [U, H| 
stress the importance of multifractality in the 
study of heart rate variability and suggest that 
it can provide a new observational window into 
the complexity mechanism of heart rate control. 
The study also highlights the need for evaluat- 
ing new nonlinear parameters for a better phys- 
iologcal investigation and for finding new clini- 
cal applications. Here we present a novel auto- 
mated scheme to compute the f(a) spectrum of 
a multifractal from its time series. We show that 
the scheme can be applied to synthetic as well 
as practical time series involving noise. It also 
provides us with an additional set of two inde- 
pendant parameters apart from the conventional 
OL m in and a max to characterise any general f(a) 
curve. The utility of these parameters from the 



point of view of diagnostic measures is also pur- 
sued by analysing a few class of physiological time 



series. 
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I. INTRODUCTION 

Multifractal sets and objects form the supporting 
structure of nonlinear phenomena, prime examples be- 
ing strange attractors of chaotic dynamical systems 0, E| 
regions of high vorticity in fully developed turbulence 
[ETE 0] and fractal growth patterns || Q . Multifractal 
analysis has also been applied in a variety of other fields, 
such as, to describe the morphologic and hydrologic char- 
acteristics of river basins [lO, [ll[ and to analyse the ve- 
locities of solar wind plasma in the inner heliosphere [l^] • 
In the case of chaotic dynamical systems, their long time 
evolution takes place in a subset of the phase space called 
strange attractors, which are characterised by a spectrum 
of generalised dimensions D„ and the associated sin- 
gularity spectrum f(a) [3, OEll- The f(a) spectrum 
provides a mathematically precise and naturally intuitive 
description of the multifractal measure in terms of inter 
woven sets with singularity strength a, whose fractal di- 
mension is f(a). 

For simple chaotic systems, such as one dimensional 
maps, f(a) spectrum can be determined analytically. To 
evaluate the /(a) spectrum from the time series, there 
are basically two methods. In the conventional method 
H3, EH, one nrs t computes the D q spectrum from the 
time series and use the fact that the transformation from 
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D q to f(a) is a Legendre transformation determined by 
the equations [H, Qjl] 

a = ^[( ? -l)i? g ] (1) 

/(a) =qa-(q- l)D q (2) 

However, such a procedure is generally considered to be 
very difficult when done subjectively, as it involves first 
smoothing the D q curve and then Legendre transforming. 
Moreover, the error bar from the smoothing procedure 
makes the estimation of f(a) more difficult and often, 
the complete spectrum cannot be recovered. 

An alternative method has been proposed in the liter- 
ature by Chhabra and Jensen [CJ] [20] for the evaluation 
of f(a) from a time series without resorting to the inter- 
mediate Legendre transform. In this method /(a) is com- 
puted directly from the slopes by plotting the normalised 
measures defined through probabilities as a function of 
logarithm of box size for different q values. The method 
gives good results and is devoid of the difficulty of Legen- 
dre transforming in the conventional method. But here 
one needs to use an optimal covering of the measure and 
the subjective evaluation of the slopes can give rise to 
error bar directly in the f(a) curve. Both these meth- 
ods generally compute only an incomplete f(a) spectrum 
(usually the top region) and the problem becomes worse 
for attractors of more than one dimension. 

Here we propose an algorithmic approach to overcome 
these difficulties and compute the complete spectrum 
from the time series for any dimension. Our scheme 
is based on the idea that the typical convex profile of 
the /(a) spectrum can be fitted by an analytic function 
involving a set of parameters. This function can be in- 
verted using Eqs. (HJ and ||2J| to get a smooth D q curve, 
which can in turn be fitted to the D q spectrum computed 
from the time series. By changing the parameters, the 
statistically best fit D q curve is chosen from which the 
final /(a) spectrum can be evaluated. It should be noted 
that we are not proposing any new method, but a new 
algorithmic approach to improve the existing methods. 

Our approach has several new features. It avoids many 
of the sources of error in the conventional method such 
as smoothing the D q curve and using the polynomial fit 
to recover the complete f(a) curve. Moreover, the whole 
procedure is made into an automated algorithmic scheme 
in the sense that once the time series is given, the scheme 
computes the D q and f(a) curves for the required embed- 
ding dimension without requiring any intermediate sub- 
jective analysis. Though the scheme is illustrated here 
for the conventional method, it can in principle be ap- 
plied for the CJ method as well. The only difference 
is that the fit has to be performed directly on the f(a) 
spectrum computed from the time series rather than the 
D q spectrum. 

Apart from the computation of the f(a) curve, another 
important outcome of our algorithmic approach is the re- 
sult that any f(a) curve can, in general, be completely 
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FIG. 1: The D q values, with error bar, of the strange attractor 
at the period doubling accumulation point of the logistic map 
calculated from the time series with 10000 data points are 
shown in the upper panel. To show the accuracy of fitting, 
these values are again plotted in the lower panel without error 
bar (dashed lines) along with the best fit curve (continuous 
line). 

characterised with the help of four independant parame- 
ters including the conventional a m i n and a max . From a 
practical point of view, this presents us with more options 
for representing the changes in the multifractal character 
of a system. Some preliminary results in this regard ob- 
tained from the analysis of a class of physiological time 
series is also included in the paper. 

Our paper is organised as follows: the algorithmic 
scheme is discussed in detail in Sec. II. The scheme is then 
tested in Sec. Ill using the time series from logistic attrac- 
tor at the period doubling accumulation point where, the 
f(a) curve is known theoretically. It is then applied to 
some other standard chaotic attractors in higher dimen- 
sions. Sec. IV considers the application of the scheme to 
practical time series, where the effect of noise on the /(a) 
spectrum is also studied. Discussions and conclusions are 
given in Sec.V. 

II. ALGORITHMIC SCHEME 
A. Computation of D q 

As the first step, the spectrum of generalised dimen- 
sions D q are computed from the time series using the 
delay embedding technique [2lj |. For this, an embedded 
space of dimension M is constructed from the scalar time 
series s(tj) as 

xl = [s(ti), s(U + r), s(ti + (M - 1)t)] (3) 
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FIG. 2: (a)The /(a) spectrum of the logistic attractor ob- 
tained directly from the D q values which is incomplete. 
(b)The f(a) spectrum computed from the best fit curve (con- 
tinuous line) in the previous figure along with the theoretical 
curve (dashed line). The agreement between the two is evi- 
dent. 

where r is a suitably chosen time delay. The generalised 
correlation sum C q (R) is given by the relative number of 
data points within a distance R from a particular (i th ) 
data point, say Ci(R), raised to the power of (q — 1) and 
averaged over N c randomly selected centres: 

N v 

C ^ = W E H(R-\x' l -x' J \) (4) 

" ■ 1 



recently proposed an algorithmic scheme which non sub- 
jectively chooses an appropriate region [22]. The scheme 
which was developed to compute D 2 , uses M-cubes in- 
stead of M-spheres and chooses only those centres where 
the M-cubes are inside the attractor region, thus avoiding 
the "edge effect" . For large R, the number of such cen- 
tres decrease and by the condition that at least iV„/100 
centres are used for the computation, a maximum value 
of R, Rmax is obtained. To avoid the region dominated 
by counting statistics only results from R > Rmin are 
taken into considerations where N V C2(R) > 10, which 
ensures that on the average at least ten data points are 
being considered. The above estimation of the scaling 
region R m in < R < Rmax is adequate for computation 
of £>2 and the results obtained matches with theoretical 
values [22]. However, for the generalised correlation di- 
mensions, there is an additional error which may occur 
for finite data sets. For large absolute values of q, the 
average over the randomly chosen centres Eq. ^ may be 
dominated by a few centres which have either large (for 
q > 0) or small (for q < 0) values of Ci(R). This biases 
the result towards a few centres which may be due to 
statistical fluctuations. This effect is particularly strong 
when q < 0, where centres with statistically small values 
of Ci(R) are the main contributors to C q (R). This is over- 
come by demanding that at least 1/10 of the centres have 
a value of Ci(R) q ~ 1 greater than C q (R). This restricts the 
range of R further and often a suitable range of R is not 
available for small values of q. We fit a straight line to 
log C q (R) versus log R for the range of R that satisfy the 
above criteria and estimate D q from the slope of the best 
fit line. The standard error on the slope is used as an 
estimate of the error on D q . 



B. Computation of /(a) 



c q {R) = w Y.^ R ) q 



(•5) 



where N v is the number of vectors. Then the spectrum 
of dimensions are given by 



D q = 



1_ Hm log C q (R) 
R^o log R 



1 



(6) 



In practical considerations, D q is computed by taking 
the slope of log C q (R) versus logi? over a region of R 
where the slope is nearly a constant i.e. a scaling re- 
gion. For large R (apart from the formal breakdown of 
the limit R — > 0), the slope does not represent D q since 
the M-spheres may extend outside the attractor, an ef- 
fect known as "edge effect" . For small R, the correlation 
sum is affected by counting statistics due to the finite 
length of the data stream. In general, the appropriate 
scaling region where both these effects are negligible is 
often chosen by subjective visual inspection. We have 



Attempting to compute the f(a) spectrum directly 
from the D q values using Eqs. (|TJ) and © leads to an in- 
complete /(a) spectrum (see Fig. [5^). This is mainly due 
to the fact that the errors in the calculation of D q makes 
the Legendre transforming numerically impractical be- 
cause of reversal of slopes. The conventional method is 
to either smoothcn the D q values or use a polynomial 
fit to recover the complete /(a) curve. Both can lead 
to large errors as has already been discussed by many 
authors. Here we follow a different procedure as given 
below. 

The f(a) function is a single valued function defined 
between the limits of a m i n and a max . Since the deriva- 
tive /'(a) = df(a)/da = q is also single valued, it fol- 
lows that f(a) has a single extremum (i.e. a maximum). 
Moreover, f(a min ) = f{a max ) = and f'(a min ) and 
f'i&max) tend to oo and — oo respectively. A simple func- 
tion which can satisfy all the above necessary conditions 
is 



f(a) = A(a 



P 1 



{a r , 



a 



,72 



(7) 



4 



Logistic Attractor 




Rossler Attractor 
a = b = 0.2,c = 7.8 



0.2 0.4 




FIG. 3: The variation of the f(a) spectrum of the logistic 
attractor with number of data points. The spectrum almost 
completely coincide for 10000 data points (dashed line) and 
5000 data points (solid line), but shows slight deviation as the 
number of data points are reduced to 3000 (thick line). 



Rossler Attractor 
a = b = 0.2,c = 7.8 



FIG. 4: The D q values (points) and its best fit curve (con- 
tinuous line) for the Rossler attractor with 10000 data points 
computed using the scheme. 



where A, 71, 72, a m i n and a max are a set of parame- 
ters characterising a particular f(a) curve. We will show 
that out of these five parameters, only four are inde- 
pendant which can unambiguously fix any general /(a) 
curve. From Eq. |2]), we get 



da 



/(«) 



(8) 



FIG. 5: The f(a) spectrum of the Rossler attractor computed 
from the best fit curve in Fig. Q. 



Substituting for f(a) from above and simplifying 



/(«) 



7i 



72 



(9) 



The form of /(a) assumed in Eq. ([7]) implies that for 
it to be a well behaved function, A should be positive 
and 71,72 > 0. If 71,72 < 0, then f(a) — > 00, as a — > 
drain , Oi ma x ■ Imposing the condition that the slope of 
f(a) (that is, j^f(a)) should be 00 at a = a min ,a max , 
we find from Eq. (O that this is possible only if both 
71)72 < 1- Thus the range of 71,72 should be restricted 
to 



< 71 , 72 < 1 



(10) 



Corresponding to q = 1, there is a value of a(= a\) and 
/(«)(= f(a±)) such that 



Di = a 1 = f(ai) 
Putting q — 1 in Eq. we get 



«i 



7i 



12 



ai 



(11) 



(12) 



Using ai, a m in, ct m ax and 71 as input parameters, 72 can 
be calculated from this equation. These values are then 
used to calculate the parameter A from the original f(a) 
fit (Eq. ([7])) with a = a\(= f(ax)). Thus only four inde- 
pcndant parameters are required to fix the f(a) curve. 

The scheme first takes ai(= Di),a min (= D^) and 
a max (= D-oo) as input parameters from the computed 
D n values and choosing an initial value for 71 in the range 
[0, 1], the parameters 72 and A are calculated. The f(a) 
spectrum is then computed in the range [a m i n ,a max ]. 



5 




Lozi ALLractor 
a=1.7,b=0.5 



FIG. 6: The upper panel shows the D q values and its best fit 
curve (continuous line) for the Lozi attractor obtained by ap- 
plying our numerical scheme. The corresponding f(a) curve 
is shown in the lower panel. Note that the D q branch for 
q > is almost flat resulting in a highly asymmetric /(a) 
curve. 



From this, the D q versus q curve is then computed by 

and fitted to the D q values 



inverting Eqs. (1} and |j2]) 
computed from the time series 



A % 2 fitting is undertaken by changing the parameter 
71 (which in turn changes 72 and A) untill the functional 
fit matches with the D q values from the time series in the 
best possible manner as indicated by the minimum of the 
X 2 value. The final f(a) spectrum is computed from the 
functional fit. 



In the scheme, D\ is used as one of the input parame- 
ters as it is obtained directly from the D q spectrum com- 
puted from the time series. But to represent the changes 
in the f(a) curve, it is more convenient to use 71 and 
72 along with a m i n and a rnax as the free parameters. 
While a m i n and a max fix the end points of the spectrum 
(at which the slope becomes 00), the values of 71 and 
72 completely specify the nature of the f(a) profile. For 
example, as the difference between 71 and 72 increases, 
the spectrum becomes more and more asymmetric be- 
tween the two branches. Moreover, they also determine 
the peak value of the spectrum D - If 71 ~ 72 ~ 0, 
then Do <C M and on the other hand, if 71 ~ 72 ~ 1, 
then Dq ~ M. Thus the four parameters a m i n , a ma x,li 
and 72 can, in principle, completely characterise any f{a) 
spectrum. 



III. APPLICATION TO SYNTHETIC DATA 

To test our scheme, it is first applied on the time series 
from the logistic attractor at the period doubling accu- 
mulation point, where the f{a) spectrum is known the- 
oretically. The analysis is done with 10000 data points. 
The D q spectrum is first computed using Eq. ^ (with 
M = 1), for q values in the range [—20, +20]. The compu- 
tation is done taking a step width of Aq = 0.1. Choosing 
-D_20i Di and Z?20 as the input values for the f{a) func- 
tion Eq. the parameters 71 and 72 are scanned in the 
range [0, 1] untill the functional fit matches the D q values 
as indicated by the x 2 minimum. Since the error in D q 
generally bulges as q — > —20, the error bar is also taken 
into account in the fitting process. The D q values with 
error bar and the best fit curve are shown in Fig. [TJ The 
complete f(a) spectrum computed from the best fit D q 
curve is shown in Fig.[2j For one dimensional maps, the 
f(a) curve can be determined theoretically To make 
a comparison, the theoretical /(a) curve is superimposed 
on the computed one in Fig. [2] Also shown in Fig. [2^ is 
the incomplete /(a) spectrum computed directly from 
the D q values. 

Since the /(a) spectrum for the logistic attractor is 
exactly known, it can also be used to test our scheme 
with respect to the number of data points required in a 
time series for a reasonable estimate of the f{a) spec- 
trum. This is shown in Fig. [31 where the spectrum for 
the logistic attractor for three different number of data 
streams, namely, 3000, 5000 and 10000 are shown. It 
turns out that, a reasonable approximation to the f(a) 
spectrum can be obtained atleast with 3000 data points 
using our scheme in one dimension. 

As the second example, we use time series generated 
from another standard chaotic attractor, namely the 
Rossler attractor, for parameter values a = 0.2, 6 = 
0.2, c = 7.8, with a time step of At = 0.1. The D q and 
f(a) spectrum are computed as above taking the total 
number of 10000 data points. The D q values and best 
fit curves for are shown in Fig. SI while the f{a) curve is 
shown in Fig. [5] 

In the above two cases, the f{a) curves are almost 
symmetric between the two branches. We now show that 
our scheme is also useful for the computation of more 
general types of /(a) curves. An example is that of the 
standard Lozi attractor. It is found that the D q spectrum 
in this case is almost flat for q > 0. But the D q spectrum 
can be accurately fitted by using the standard form of the 
f(a) function, as shown in Fig. The corresponding 
f(a) spectrum computed from the best fit D q curve is 
also shown in Fig. [6] (lower panel), which turns out to 
be highly asymmetric between the two branches. While 
the parameter values 71 and 72 for the Rossler attractor 
are very close (71 = 0.65,72 = 0.60) as reflected in the 
nearly symmetric f(a) profile, that for Lozi attractor are 
completely different with 71 = 0.09 and 72 = 0.60. The 
scheme has also been applied to compute the spectrum 
of other standard chaotic attractors, such as, Henon and 
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FIG. 7: The D q spectrum for white noise and red noise from 
10000 data points. Note that the latter behaves like a chaotic 
system with a well defined D q curve. 
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FIG. 8: The /(a) spectrum for white noise and red noise. 
The scheme accurately determines the /(a) spectrum of white 
noise which is expected to be a 8 function corresponding to 
the embedding dimension. 

Lorenz. 



IV. APPLICATION TO REAL WORLD DATA 

Before the scheme is applied to practical time series, 
it is important to test it with time series involving noise. 
Since pure white noise is scale free, one expects the cor- 



FIG. 9: Result of addition of white noise to the data from 
the logistic attractor. The f(a) spectrum for a noise level 
of 20% (dotted line) and 50% (solid line) are shown along 
with that of logistic attractor (points). The number of data 
points used for computation are 10000. Note that for 50% of 
noise contamination, the peak of the spectrum (Do) almost 
touches the embedding dimension 1. For comparison, the 
f(a) spectrum of pure random noise in one dimension is also 
shown (dashed line). 

responding /(a) spectrum to be ideally a S function with 
/(a) = a — M, the embedding dimension. In Fig. we 
show the D q spectrum of pure white noise and a colored 
noise with spectral index 2.0 (called red noise), computed 
from the respective time series for embedding dimension 
M = 3. The /(a) spectrum computed for both are shown 
in Fig. [51 While the f(a) spectrum of white noise is a 
S function as expected, that of colored noise looks like a 
normal f(a) spectrum. 

To study the effect of noise on the f(a) spectrum of a 
chaotic attractor, we generate two time series by adding 
20% and 50% white noise to the data from the logistic 
attractor using 10000 data points. White noise is used 
since it is scale free and hence can significantly alter the 
/(a) spectrum, while colored noise behaves much like a 
chaotic attractor with a well defined spectrum. Fig. [9] 
shows the /(a) spectrum of the logistic attractor added 
with 20% and 50% white noise, along with that of the lo- 
gistic attractor. As the percentage of noise increases, the 
spectrum tends more and more towards a delta functon, 
centered around the embedding dimension M = 1. To 
get a proper comparison, the figure also shows the f(a) 
spectrum of pure random noise in one dimension. 

We now apply our scheme to some real world data. We 
choose time series from two important fields, namely, as- 
trophysics and physiology, where methods and concepts 
from nonlinear dynamics are constantly being applied. 
The first example is the X-ray light curve from a promi- 
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FIG. 10: A part of the light curves from two temporal states 
of the black hole system GRS1915+105. 



FIG. 12: The f(a) spectrum for the black hole states com- 
puted from the best fit D q curves in the previous figure. 
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FIG. 11: The D q spectrum along with the best fit curve for 
the two states 6 and k of the black hole system for embedding 
dimension M = 3. 



FIG. 13: A part of the EEG and ECG signals analysed in this 
work. 



nent back hole binary, GRS1915+105. The light curves 
from this black hole system have been classified into 12 
temporal states by Belloni et.al [23| based on RXTE ob- 
servation. Here we choose data from two representative 
classes, 9 and n, and generate continuous light curves of 
approximately 6000 data points for both class. Fig. [TU1 
shows a part of the light curves used for the analysis. 
Using surrogate analysis, we have recently shown [HI HE] 
that light curves from more than half of the 12 tempo- 
ral states (including 9 and k) show significant deviation 



from stochastic behavior. The saturated value of cor- 
relation dimension for both are < 3. Hence M — 3 is 
chosen for applying our numerical scheme. The result of 
applying our scheme to the two light curves are shown 
in Fig. [11] and Fig. [T2l The former shows the computed 
D q spectrum along with the best fit curve, while the lat- 
ter shows the /(a) spectrum computed from the best fit 
curves. Note that, though the D q values are discontinu- 
ous for large negative q values, one can statistically fit a 
smooth curve for D q from which the f(a) spectrum can 
be derived. 
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FIG. 14: The variation of correlation dimension Di (with 
error bar) as a function of M for the EEG and ECG signals. 



FIG. 15: The D q spectrum and the best fit curves for the 
EEG and ECG signals for M = 3. 



As the second example, we use two data sets 
from physiology, namely an EEG data and an ECG 
data. The EEG data was downloaded from the 
website of the Department of Epileptology, Univer- 
sity of Bonn while the ECG data was obtained from 
|http:/ /www.physionet.o rg/physi obank/archives] The 
EEG data is from an epileptic patient during seizure ac- 
tivity. The data consists of continuous data streams of 
about 24 seconds long and consisting of approximately 
5000 data points. The ECG data was recorded from a 
heart patient with a congestive disorder and consists of 
continuous data streams of 5400 data points with a sam- 
pling time of 0.04 seconds. Both signals are shown in 
Fig. [T3J First we compute the correlation dimension D2 
of both signals by applying the nonsubjective scheme [22j 
recently proposed by us, with the result shown in Fig.[T4l 
For both signals, D2 saturate well below M — 3. The re- 
sults of applying our f(a) scheme are shown in Fig. [15] 
and Fig. 1 1 61 with multifractal character evident in both 
cases. Thus it is clear that the scheme can be success- 
fully employed to compute the D q and /(a) spectrum 
from practical time scries of finite data streams even with 
noise contamination, provided there exists an underlying 
chaotic attractor. 

Finally, we would also like to stress the importance 
of computing the f{a) spectrum using an automated 
scheme, such as the one presented here. In the four ex- 
amples of real world data that we have analysed above, 
the saturated D2 values are approximately same, but the 
f(a) spectra are quite different. The subtle changes in 
the /(a) spectra can be better studied using the present 
scheme. This is because, the scheme provides an ad- 
ditional set of two independant parameters 71 and 72, 
apart from a m i n and a max which typically characterise 




FIG. 16: The /(a) spectrum for the EEG and ECG signals 
computed from the best fit D q curves in the previous figure. 



the changes in the f(a) profile. The utility of these pa- 
rameters can be seen from Table U which shows the 
result of a preliminary study made on some limited phys- 
iological time series. We have analysed four class of phys- 
iological data sets downloaded from the above mentioned 
websites. They are EEG and ECG signals from healthy 
human beings, EEG signal during epileptic seizure and 
ECG signals from patients with congestive heart failure. 
We have analysed five data streams for each class and 
Table [I] shows the average values of the two parame- 
ters 71 and 72 for each class along with a m i n and a rnax , 
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TABLE I: The parameter values computed by our scheme for physiological data sets corresponding to four different class. The 
average values of five data streams for each class are shown. 

Data Class a min a max 71 72 I71 - 72 1 

EEG 

Healthy 1.71 ± 0.08 4.18 ±0.17 0.37 ±0.10 0.32 ±0.14 0.05 

EEG 

Epileptic Seizure 1.28 ± 0.06 2.85 ±0.14 0.34 ± 0.08 0.05 ± 0.03 0.29 

ECG 

Healthy 1.46 ±0.12 4.30 ± 0.22 0.69 ± 0.08 0.58 ± 0.05 0.11 

ECG 

Congestive Heart Failure 1.88 ±0.10 4.07 ±0.18 0.29 ± 0.06 0.08 ± 0.04 0.21 



with error bar showing the range of variation. We find 
that the values of I71 — 72 1 are significantly different for 
the healthy data and those with physiological disorder in 
both cases. Thus, just like \(x max — a m i n \ which is con- 
ventionally used to characterise the geometric complexity 
of a multifractal, the quantity \ ji — 72 1 may also be a po- 
tential candidate to quantify the inherent changes in the 
multifractal character of the system. Ofcourse, these re- 
sults are only preliminary and have to be confirmed using 
a much larger number of data sets and requires an exten- 
sive analysis. Nevertheless, the initial results indicate 
that our scheme can be efficiently employed to quantify 
the changes in the multifractal character between various 
states of a complex system (such as the black hole sys- 
tem) or track changes in the f(a) spectrum arising out 
of various physiological disorder as reflected in the ECG 
or EEG time series, as shown here. We hope the scheme 
can be better utilised in this regard. 



V. DISCUSSION AND CONCLUSION 

Computing the multifractal spectrum of a chaotic at- 
tractor from its time series is generally considered to be 
a difficult task. Often, only a part of the f(a) spectrum 
can be recovered numerically from the D q curve due to 
various reasons such as, the errors in the computation 
of D q and the Legendre transforming involved. Here we 
show that the existing methods can be improved using 
a new algorithmic approach by which the complete f(a) 
spectrum can be evaluated from a time series. 

The scheme first assumes an analytical function for 
the /(a) curve, from which a functional fit for the com- 
puted D q values can be obtained by the inverse Legendre 
transform. The best fit curve is then used to derive the 
complete f(a) spectrum. The scheme is illustrated using 
time series from standard low dimensional chaotic sys- 
tems and then applied to a variety of practical data from 
real world. 



We have recently shown [26| that the the /(a) gives 
information only upto two scales evenif the underlying 
multiplicative process involves more than two. As a con- 
sequence, the f(a) spectrum of a chaotic attractor can, 
in general, be mapped onto that of a two scale Cantor 
set. This also provides an alternative method for com- 
puting the f(a) spectrum and to characterise a chaotic 
attractor in terms of the three independant parameters 
of a two scale Cantor set. 

In contrast, an important aspect of the present scheme 
is that, an analytic function is proposed (which is proba- 
bly unique) to fit all convex f(a) curves in general. Since 
the whole process is automated and the analysis is done 
under identical conditions prescribed by the algorithmic 
scheme, the resulting parameters characterising the spec- 
trum (derived from a fitting function for f(a)) give a 
better representation for comparison between data sets. 
This is especially important in the case of real world data 
since the changes in the same system, such as for exam- 
ple, due to some changes in the parameter, can be com- 
pared in a non-subjective manner. 

In order to show this explicitely, the analysis of a few 
class of physiological data is also presented. We find 
that the subtle variations in the f(a) curve can be bet- 
ter characterised using our algorithmic scheme. Though 
our results indicate that the parameters may also be use- 
ful from a diagnostic point of view, this requires a much 
more comprehensive analysis using large number of data 
sets for confirmation. This is currently under way and 
will be presented elsewhere. 
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